Shape Variables and Shape Spaces
Antigoni Kaliontzopoulou, CIBIO/InBIO, University of Porto
15 October, 2019
From GPA to Shape Data

- In this lecture we:
- Obtain variables we can use for analyses
- Review the relative position of individuals in shape space
- Visualize shape variation
Shape Spaces from GPA
- After GPA, each landmark configuration is a point in shape space
- The ācenterā of this universe is the consensus (global mean shape)
- Each shape occupies a unique point in shape space
- The metric of shape space is the Procrustes distance \(\small{D}_{Proc}\)
\[\tiny{D}_{Proc}=\sqrt{\sum_{i,j}^{k,p}\left(\mathbf{Y}_{1.ij}-\mathbf{Y}_{2.ij}\right)^2}\]
Shape Spaces from GPA
- After GPA, each landmark configuration is a point in shape space
- The ācenterā of this universe is the consensus (global mean shape)
- Each shape occupies a unique point in shape space
- The metric of shape space is the Procrustes distance \(\small{D}_{Proc}\)
\[\tiny{D}_{Proc}=\sqrt{\sum_{i,j}^{k,p}\left(\mathbf{Y}_{1.ij}-\mathbf{Y}_{2.ij}\right)^2}\]

Kendallās Tangent Space Coordinates
- Orthogonal projection of shapes to linear tangent space
\[\small\mathbf{X}'=\mathbf{X\left(I_{kp}-X_c^T(X_cX_c^T)^{-1}X_c\right)}\]

- Specimen scores from projection are Kendallās tangent space coordinates (Dryden and Mardia, 1993; Kent 1994)
- Perform PCA and retain only 2p-4 (3D: 3p-7) shape variables with variation if desired (NOTE: permutation via RRPP is now standard practice, so this step is not strictly necessary!)
Kendallās Tangent Space Coordinates
- Random uniform shapes ARE uniformally distributed in Tangent Space!

Same set of 2,000 randomly generated triangles
The Thin-Plate Spline
- TPS is a smooth interpolation function (i.e., a doubly-differentiated equation)
- It models the differences in landmark locations between the reference and target specimen for each coordinate dimension separately (x,y)

The Thin-Plate Spline
- TPS is a smooth interpolation function (i.e., a doubly-differentiated equation)
- It models the differences in landmark locations between the reference and target specimen for each coordinate dimension separately (x,y)

General TPS Model for \(\small{1}\) landmark: \[\tiny\begin{bmatrix} x^1\\y^1\end{bmatrix}=A\begin{bmatrix}1\\x\\y\end{bmatrix}+\sum_{i=1}^pw_iU(r_i)\]
where \(\tiny{A}\begin{bmatrix}1\\x\\y\end{bmatrix}\) represents the affine (uniform) component and \(\small\sum_{i=1}^pw_iU(r_i)\) represents to non-affine (non-uniform) component
TPS Computations
- 1: From the reference (consensus) estimate: \(\tiny\mathbf{L}= \left[ \begin{array}{c|c} \mathbf{P} & \mathbf{Q}\\ \hline \mathbf{Q}^T & \mathbf{0} \\ \end{array} \right]\)
where: \(\tiny\mathbf{Q}= \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ 1 & x_p & y_p \end{bmatrix}\) & \(\tiny\mathbf{0}= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\)
\(\tiny\mathbf{P}=\tiny\begin{bmatrix} 0 & U(r_{12}) & U(r_{12}) & \dots & U(r_{1p}) \\ U(r_{21}) & 0 & U(r_{23}) & \dots & U(r_{2p}) \\ U(r_{31}) & U(r_{32}) & 0 & \dots & U(r_{3p}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ U(r_{p1}) & U(r_{p2}) & U(r_{p3}) & \dots & 0 \end{bmatrix}\)
TPS Computations
- 1: From the reference (consensus) estimate: \(\tiny\mathbf{L}= \left[ \begin{array}{c|c} \mathbf{P} & \mathbf{Q}\\ \hline \mathbf{Q}^T & \mathbf{0} \\ \end{array} \right]\)
where: \(\tiny\mathbf{Q}= \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ 1 & x_p & y_p \end{bmatrix}\) & \(\tiny\mathbf{0}= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\)
\(\tiny\mathbf{P}=\tiny\begin{bmatrix} 0 & U(r_{12}) & U(r_{12}) & \dots & U(r_{1p}) \\ U(r_{21}) & 0 & U(r_{23}) & \dots & U(r_{2p}) \\ U(r_{31}) & U(r_{32}) & 0 & \dots & U(r_{3p}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ U(r_{p1}) & U(r_{p2}) & U(r_{p3}) & \dots & 0 \end{bmatrix}\)
Here, \(\small\mathbf{Q}\) contains the coordinates of the reference and \(\small\mathbf{P}\) is found as: \(\small{U(r_{ij})}=r_{ij}^2ln(r_{ij})\) with \(\small{r_{ij}}\) as the distance between the \(\small{i^{th}}\) and \(\small{j^{th}}\) landmarks
NOTES: \(\small{0}\) is \(\small{4 \times 4}\) for 3D data. Also, for 3D data: \(\small{U(r_{ij})}=|{r_{ij}}|\)
TPS Computations (Cont.)
- 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
TPS Computations (Cont.)
- 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
- 3: Perform eigenanalysis of \(\small\mathbf{L}_p^{-1}\) to obtain Principal Warps which are the \(\small{p-3}\) eigenvectors containing variation (\(\small{p-4}\) for 3D)
\[\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\]
TPS Computations (Cont.)
- 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
- 3: Perform eigenanalysis of \(\small\mathbf{L}_p^{-1}\) to obtain Principal Warps which are the \(\small{p-3}\) eigenvectors containing variation (\(\small{p-4}\) for 3D)
\[\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\]
- 4: Estimate non-uniform shape variables (Partial Warp Scores) via projection onto principal warps \(\small\mathbf{E}\):
\[\small\mathbf{W}=\mathbf{V(I_2\otimes{E})=V}\begin{bmatrix}\mathbf{E} & \mathbf{0} \\ \mathbf{0} & \mathbf{E} \end{bmatrix})\]
where \(\small\mathbf{V}=\mathbf{[V_x|V_y]}\) with \(\small\mathbf{V_x}=\mathbf{(X-\overline{X}_{ref})}\) & \(\small\mathbf{V_y}=\mathbf{(Y-\overline{Y}_{ref})}\)
TPS Computations (Cont.)
- 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
- 3: Perform eigenanalysis of \(\small\mathbf{L}_p^{-1}\) to obtain Principal Warps which are the \(\small{p-3}\) eigenvectors containing variation (\(\small{p-4}\) for 3D)
\[\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\]
- 4: Estimate non-uniform shape variables (Partial Warp Scores) via projection onto principal warps \(\small\mathbf{E}\):
\[\small\mathbf{W}=\mathbf{V(I_2\otimes{E})=V}\begin{bmatrix}\mathbf{E} & \mathbf{0} \\ \mathbf{0} & \mathbf{E} \end{bmatrix})\]
where \(\small\mathbf{V}=\mathbf{[V_x|V_y]}\) with \(\small\mathbf{V_x}=\mathbf{(X-\overline{X}_{ref})}\) & \(\small\mathbf{V_y}=\mathbf{(Y-\overline{Y}_{ref})}\)
Note, one can weight partial warp scores inversely by spatial scale as: \(\small\mathbf{W}=\mathbf{V(I_2\otimes{E\Lambda^{-\alpha/2}})}\)
- There are \(\small{2p-6}\) partial warp scores (\(\small{3p-12}\) for 3D data)
Total Shape Variation
- Shape variation includes BOTH uniform and non-uniform components

Together, the Partial Warps + Uniform component scores comprise the total set of shape variables: \(\small\mathbf{W}=\mathbf{V(I_2\otimes{E})|U}\)
There are \(\small{2p-4}\) shape variables for 2D data and \(\small{3p-7}\) for 3D data
TPS Model: Example
- Model shape differences from a square to a kite (only non-uniform shape differences present)

\(\tiny\mathbf{X}_{sqr}=\begin{bmatrix} 0 & 0.5\\ -0.5 & 0\\ 0 & -0.5\\ 0.5 & 0 \end{bmatrix}\) and \(\tiny\mathbf{X}_{kite}=\begin{bmatrix} 0 & 0.375\\ -0.5 & 0.125\\ 0 & -0.625\\ 0.5 & 0.125 \end{bmatrix}\)
TPS Model: Example (Cont.)
- Obtain: \(\tiny\mathbf{L}= \left[ \begin{array}{c|c} \mathbf{P} & \mathbf{Q}\\ \hline \mathbf{Q}^T & \mathbf{0} \\ \end{array} \right]=\left[ \begin{array}{cccc|ccc} 0 & -0.173 & 0 & -0.173 & 1 & 0 & 0.5\\ -0.173 & 0 & -0.173 & 0 & 1 & -0.5 & 0 \\ 0 & -0.173 & 0 & -0.173 & 1 & 0 & -0.5\\ -0.173 & 0 & -0.173 & 0 & 1 & -0.5 & 0 \\ \hline 1 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & -0.5 & 0 & 0.5 & 0 & 0 & 0 \\ 0.5 & 0 & -0.5 & 0 & 0 & 0 & 0 \\ \end{array} \right]\)
- Then calculate: \(\tiny\mathbf{L^{-1}}=\left[ \begin{array}{cccc|ccc} 0.7271 & -0.7271 & 0.7271 & -0.7271 & 0.25 & 0 & 1\\ -0.7271 & 0.7271 & -0.7271 & 0.7271 & 0.25 & -1 & 0 \\ 0.7271 & -0.7271 & 0.7271 & -0.7271 & 0.25 & 0 & -1\\ -0.7271 & 0.7271 & -0.7271 & 0.7271 & 0.25 & 1 & 0 \\ \hline 0.25 & 0.25 & 0.25 & 0.25 & 0.087 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 & 0 \\ \end{array} \right]\)
- \(\tiny\mathbf{L_p^{-1}}=\left[ \begin{array} 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ \end{array} \right]\)
TPS Model: Example (Cont.)
- Solve: \(\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\)
\(\tiny\mathbf{L_p^{-1}}=\left[ \begin{array} 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ \end{array} \right]=\)
- \(\tiny\mathbf{E}=\left[ \begin{array} \mathbf{0.5} & 0.866 & 0 & 0 \\ \mathbf{-0.5} & 0.289 & 0 & -0.0816 \\ \mathbf{0.5} & -0.289 & 0.707 & -0.408 \\ \mathbf{-0.5} & 0.289 & -0.707 & 0.408 \\ \end{array} \right]\) & \(\tiny\mathbf{\Lambda}=\left[ \begin{array} \mathbf{2.885} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array} \right]\)
- Obtain shape variables: \(\small\mathbf{W=VE}\) using: \(\tiny\mathbf{V_{kite}}=\begin{bmatrix} 0 & 0 & 0 & 0 & -1.25 & -1.25 & -1.25 & -1.25 \end{bmatrix}\)
- \(\tiny\mathbf{W=VE}=\mathbf{V_{kite}}\begin{bmatrix} 0.5 & 0 \\ -0.5 & 0 \\ 0.5 & 0 \\ -0.5 & 0 \\ 0 & 0.5 \\ 0 & -0.5 \\ 0 & 0.5 \\ 0 & -0.5 \\ \end{bmatrix}=\begin{bmatrix} 0 \\ -0.25 \end{bmatrix}\)
- NOTE: Uniform shape variables, \(\small\mathbf{U_1}=0\) & \(\small\mathbf{U_2} = 0\) for this example
Which Shape Variables to Use?
- We now have 3 shape variable options:
- GPA-aligned coordinates
- Kendallās Tangent Space coordinates (GPA + orthogonal projection)
- TPS variables (partial warp scores + uniform shape)
- Which ones should we use?
Which Shape Variables to Use?
- We now have 3 shape variable options:
- GPA-aligned coordinates
- Kendallās Tangent Space coordinates (GPA + orthogonal projection)
- TPS variables (partial warp scores + uniform shape)
- Which ones should we use?

- Plot of shape distances for Podarcis data obtained from GPA-shape space and Kendallās Tangent space. Correlation of 1.0 in this case (usually very, very high)
- AND correlation between PWS+U vs.Ā Tangent space is ALWAYS 1.0 (see Rohlf 1999)
Shape Variables: Conclusions
- Proc. Tangent coordinates & TPS shape variables yield identical results
- They are simply rotations of one another
- GPA-aligned coordinates (unprojected) are close, but not exact
- Difference due to curvature of shape space
- Use Kendall“s Tangent Space Coordinates to represent shape
- These are found from GPA+Orthogonal Projection
- NOTE: Most current papers use tangent space coordinates, aka Procrustes residuals combined with resampling tests (to accomodate singular dimensions)
- TPS then, is frequently used to visualize shape differences
Exploring Shape Variation
- GM data are points in high-dimensional shape space
- (e.g.Ā for \(\small{p=20}\) & \(\small{k=3}\), shape space comprises 53 dimensions)
- Our human brains are limited to perception in up to 3 dimensions
- We need tools for visualizing and exploring shape variation in high-dimensional data spaces
- Principal Components Analysis (PCA) is one such tool
Principal Components Analysis (PCA)
- PCA combines orthogonal rotation and projection to arrive at a summary plot of the dataspace
- The objective: to summarize most of the variation in as few dimensions as possible
- It consists of a rigid rotation of the data based on directions of variation, followed by projection to those new summary axes
- First one finds the set of axes the describe progressively less variation in the data
- Next, the data are rotated so that the main axis of variation (PC1) is horizontal
- Subsequent axes are orthogonal to PC1
PCA: Conceptual Visualization

PCA: Conceptual Visualization

- PCA has performed a rigid rotation of the original data, so that variation is aligned with the new (PC) axes
PCA: Standard Computations
Principal component analysis (PCA) is two things: (1) a singular-value decomposition (SVD) of a symmetric matrix and (2) projection of mean-centered or standardized data onto eigenvectors from SVD.
Using mean-centered data: \(\small\mathbf{Y}_c\) we calculate:
\[\small\hat{\mathbf{\Sigma}}=\frac{\mathbf{Y}^{T}_c\mathbf{Y}_c}{n-1}\]
We then decompose the covariance matrix via eigen-analysis (SVD):
\[\small\hat{\mathbf{\Sigma}}=\mathbf{U} \mathbf{\Lambda} \mathbf{U}^T\]
This step identifies a set of orthogonal axes describing directions of variation. These are the columns of \(\small\mathbf{U}\). Next we project the data onto these axes as:
\[\small\mathbf{P} = \mathbf{Y}_c\mathbf{U}\]
Here, \(\small\mathbf{P}\) represent the set of projection scores (principal component scores) on each of the PC axes. These are used in the subsequent plot.
Finally, the percent variation explained by any principal component (column of \(\mathbf{U}\)) is \[\frac{\lambda_i}{tr\mathbf{\Lambda}}\]
PCA: Alternative Computations
Note that one can alternatively perform singular-value decomposition (SVD) directly on \(\small\mathbf{Y}_c\):
\[\small{\mathbf{Y}_c}=\mathbf{V} \mathbf{D} \mathbf{U}^T\]
Here, \(\small\mathbf{D^2}\) expresses the percent-variation explained by each PC-axis, which are found as the right-singular vectors in \(\small\mathbf{U}\).
PC scores are found as:
\[\small\mathbf{P} = \mathbf{VD}\]
NOTE: eigenanalysis of the inner product \(\small\mathbf{Y}^{T}_c\mathbf{Y}_c\) yields \(\small\mathbf{U}\), while eigenanalysis of the outer product \(\small\mathbf{Y}_c\mathbf{Y}^{T}_c\) yields \(\small\mathbf{V}\)
PCA of Shape Data: Relative Warps
\[\small{\mathbf{W^*}}=\mathbf{V} \mathbf{D} \mathbf{U}^T\]
- Here, \(\small\mathbf{VD}\) are the PC scores (relative warp scores) of the shape data
PCA of Shape Data: Relative Warps
\[\small{\mathbf{W^{*}}}=\mathbf{V} \mathbf{D} \mathbf{U}^T\]
Here, \(\small\mathbf{VD}\) are the PC scores (relative warp scores) of the shape data
NOTE: if W found as: \(\small\mathbf{W}=\mathbf{V(I_2\otimes{E\Lambda^{-\alpha/2}})}\) PCA=RWA is weighted by spatial scale of deformations
- In this case
- \(\small\alpha = 0\): unweighted analysis
- \(\small\alpha < 0\): emphasizes small-scale variation
- \(\small\alpha > 0\): emphasizes large-scale variation
- Careful! \(\small\alpha\neq{0}\) distorts data space, so it changes the relationships between observations!
RWA: Example
PCA of Podarcis data


RWA: Example (Cont.)
- PCA using different \(\small\alpha\) values (\(\small{1}\) & \(\small{-1}\))

RWA: Thoughts
Ordination through PCA very useful approach for exploring shape space
\(\small\alpha\neq{0}\) in RWA allows emphasizing small- or large-scale deformations
However, \(\small\alpha\neq{0}\) changes shape space, so we need to choose carefully
In practice, there is generally no conceptual reason of biological relevance for emphasizing deformations at different spatial scales
\(\small\alpha{=0}\) is preferred: no a priori, arbitrary emphasis on specific variables, and it allows the use of total shape space (i.e.Ā PW scores + U)
Visualizing Shape Differences
Visualizing Shape Differences
- TPS allows us to visualize what exactly changes in shape
As an interpolation function, TPS allows us to precisely map all points of one shape to those of another shape, and visualize how they differ
- One may wish to visualize:
- Differences between groups
- Extremes of PC axes
- How shape changes with size (allometry)
- Shape evolution (differences between ancestor and descendent on a phylogeny)
One can also magnify shape differences to aid biological interpretation

*TPS of Ref->Target is a vector in shape space, which can be extended in length
Visualizing Shape Differences: Example
- TPS of male vs.Ā female Podarcis (5X magnification)


- Major difference in posterior region of skull
Shape Variables: Conclusions
- Produce shape variables through projection to tangent space:
- Kendall“s tangent space coordinates (Procrustes residuals [GPA+Projection])
- Explore shape variation through ordination in tangent space:
- PCA of shape data (aka Relative Warps)
- Visualize shape variation to make biological interpretations:
- Use TPS to produce deformation grids